/******************************************
 
Results: HD-Mixed Line vs LD-Mixed Line
 
*******************************************/

*--------------------------------------------------
* Setup: Memory, Schemes, and Working Directory
*--------------------------------------------------

* Disable pagination
set more off

* Clear all existing data and settings
clear all

* Configure memory segments and plotting scheme
set segmentsize 3g
set scheme plotplainblind


*--------------------------------------------------
* Switches:  Sections On/Off
*--------------------------------------------------
local main_tables         1
local dynamics            1
local dynamics_binsreg    1
local ritest              1

*--------------------------------------------------
* Main Regressions
*--------------------------------------------------
if `main_tables' == 1 {

    * Load Cleaned Data
    use "$Data/Final/linelevel_data.dta", clear
    sort line Date 

    /*--------------------------------------------------
    * Main Table (Line-level Output) - Table 2
    *--------------------------------------------------*/
    areg ln_output HD_Mixed i.line i.Shift, absorb(doy) cluster(line_team)
	eststo A1
    local obs1 = e(N)
    su ln_output if e(sample) == 1
    local meanA = string(round(r(mean), 0.01))
    boottest HD_Mixed, boottype(wild)
	local lb1 = r(CI)[1,1]
	local ub1 = r(CI)[1,2]
	local LBboot1 = string(round(real("`lb1'"), 0.001))
    local UBboot1 = string(round(real("`ub1'"), 0.001))
    *eststo A1
    *local obs1 = e(N)
    *su ln_output if e(sample) == 1
    *local meanA = string(round(r(mean), 0.01))
     
    areg ln_productioncfc HD_Mixed i.line i.Shift, absorb(doy) cluster(line_team)
    boottest HD_Mixed, boottype(wild)
	local lb2 = r(CI)[1,1]
	local ub2 = r(CI)[1,2]
	local LBboot2 = string(round(real("`lb2'"), 0.001))
    local UBboot2 = string(round(real("`ub2'"), 0.001))
    eststo A2
    local obs2 = e(N)
    su productionincfc if e(sample) == 1
    local meanB = string(round(r(mean), 0.01))
     
    * Export Results to LaTeX
    estout A1 A2 using "$Output/Tables/tables_main", style(tex) replace ///
        keep(HD_Mixed) ///
        cells(b(star fmt(%9.4f)) se(par)) ///
        nolabel collabels(none) mlabels(none) starlevels(* 0.10 ** 0.05 *** 0.01) ///
        varlabels(HD_Mixed "HD-Mixed Line (LD Non-Mixed)")
    
    * Create LaTeX Footnotes
	local tex "Bootstrap(Wild Cluster) C.I. & [`LBboot1', `UBboot1'] & [`LBboot2', `UBboot2'] \\ "
    local tex " `tex' \hline "
    local tex "`tex' Observations & `obs1' & `obs2'  \\"
    local tex "`tex' Mean Dep. Var & `meanA' & `meanB'  \\"
    local tex "`tex' Day Effects & Yes & Yes   \\"
    local tex "`tex' Shift Effects & Yes & Yes \\"
    local tex "`tex' Line Effects & Yes & Yes  \\"
    local tex "`tex' \multicolumn{3}{p{8cm}}{\tiny \textit{Notes:} Standard errors clustered at production line-cohort-level. *** p<0.01, ** p<0.05, * p<0.1.} \\ \end{tabular} }"
    	
    * Append Footnotes to the LaTeX Table
    esttab A1 A2 using "$Output/Tables/tables_maina", style(tex) replace booktabs ///
        d(*) nolabel collabels(none) noobs postfoot("`tex'") nonum ///
        mtitles("(1)" "(2)") ///
        mgroups("Log Output (Pieces)" "Total Output (Boxes)", pattern(1 1) ///
        prefix(\multicolumn{@span}{c}{) suffix(}) span erepeat(\cmidrule(lr){@span}))
    	
    /*--------------------------------------------------
    * Main Table B.3 (Line X Variety Fixed Effects)
    *--------------------------------------------------*/
    qui areg ln_output HD_Mixed lv* i.Shift, absorb(doy) cluster(line_team)
    boottest HD_Mixed, boottype(wild)
    eststo B1
    local obsB1 = e(N)
    su ln_output if e(sample) == 1
    local meanB1 = string(round(r(mean), 0.01))
     
    qui areg ln_productioncfc HD_Mixed lv* i.Shift, absorb(doy) cluster(line_team)
    boottest HD_Mixed, boottype(wild)
    eststo B2
    local obsB2 = e(N)
    su ln_productioncfc if e(sample) == 1
    local meanB2 = string(round(r(mean), 0.01))
     
    estout B1 B2 using "$Output/Tables/tables_main_lv", style(tex) replace ///
        keep(HD_Mixed) ///
        cells(b(star fmt(%9.4f)) se(par)) ///
        nolabel collabels(none) mlabels(none) starlevels(* 0.10 ** 0.05 *** 0.01) ///
        varlabels(HD_Mixed "HD-Mixed Line (LD Non-Mixed)")
    
    * Create LaTeX Footnotes for Table B.3
    local tex " \\ \hline"
    local tex "`tex' Observations & `obsB1' & `obsB2'  \\"
    local tex "`tex' Mean Dep. Var & `meanB1' & `meanB2'  \\"
    local tex "`tex' Day Effects & Yes & Yes   \\"
    local tex "`tex' Shift Effects & Yes & Yes \\"
    local tex "`tex' Line X Variety Effects & Yes & Yes  \\"
    local tex "`tex' \multicolumn{3}{p{8cm}}{\tiny \textit{Notes:} Standard errors clustered at production line-cohort-level. *** p<0.01, ** p<0.05, * p<0.1.} \\ \end{tabular} }"
    	
    * Append Footnotes to the LaTeX Table B.3
    esttab B1 B2 using "$Output/Tables/tables_main_lva", style(tex) replace booktabs ///
        d(*) nolabel collabels(none) noobs postfoot("`tex'") nonum ///
        mtitles("(1)" "(2)") ///
        mgroups("Log Output (Pieces)" "Total Output (Boxes)", pattern(1 1) ///
        prefix(\multicolumn{@span}{c}{) suffix(}) span erepeat(\cmidrule(lr){@span}))
    	
    /*--------------------------------------------------
    * Main Table B.4 (Line X Day Fixed Effects)
    *--------------------------------------------------*/
    qui areg ln_output HD_Mixed i.Shift, absorb(line_day) cluster(line_team)
    boottest HD_Mixed, boottype(wild)
    eststo C1
    local obsC1 = e(N)
    su ln_output if e(sample) == 1
    local meanC1 = string(round(r(mean), 0.01))
     
    qui areg ln_productioncfc HD_Mixed i.Shift, absorb(line_day) cluster(line_team)
    boottest HD_Mixed, boottype(wild)
    eststo C2
    local obsC2 = e(N)
    su ln_productioncfc if e(sample) == 1
    local meanC2 = string(round(r(mean), 0.01))
     
    estout C1 C2 using "$Output/Tables/tables_main_ld", style(tex) replace ///
        keep(HD_Mixed) ///
        cells(b(star fmt(%9.4f)) se(par)) ///
        nolabel collabels(none) mlabels(none) starlevels(* 0.10 ** 0.05 *** 0.01) ///
        varlabels(HD_Mixed "HD-Mixed Line (LD Non-Mixed)")
    
    * Create LaTeX Footnotes for Table B.4
    local tex " \\ \hline"
    local tex "`tex' Observations & `obsC1' & `obsC2'  \\"
    local tex "`tex' Mean Dep. Var & `meanC1' & `meanC2'  \\"
    local tex "`tex' Shift Effects & Yes & Yes  \\"
    local tex "`tex' Line X Day Effects & Yes & Yes  \\"
    local tex "`tex' \multicolumn{3}{p{8cm}}{\tiny \textit{Notes:} Standard errors clustered at production line-cohort-level. *** p<0.01, ** p<0.05, * p<0.1.} \\ \end{tabular} }"
    	
    * Append Footnotes to the LaTeX Table B.4
    esttab C1 C2 using "$Output/Tables/tables_main_lda", style(tex) replace booktabs ///
        d(*) nolabel collabels(none) noobs postfoot("`tex'") nonum ///
        mtitles("(1)" "(2)") ///
        mgroups("Log Output (Pieces)" "Total Output (Boxes)", pattern(1 1) ///
        prefix(\multicolumn{@span}{c}{) suffix(}) span erepeat(\cmidrule(lr){@span}))
}

*--------------------------------------------------
* Day-by-Day Dynamics 
*--------------------------------------------------
if `dynamics' == 1 {

    * Load Cleaned Data
    use "$Data/Final/linelevel_data.dta",clear

    sort line Date 
    gen hdmixd = HD_Mixed
    egen day_of_intervention = cut(int_days), at(0(20)140)
    replace day_of_intervention = 100 if day_of_intervention == 120 
	
* Create variables needed
forv i=0(20)100 {
	g ldmixed`i' = HD_Mixed==0 & day_of_intervention==`i'
	g hdmixd`i' = HD_Mixed==1 & day_of_intervention==`i'
}

* Day-by-day dynamics
foreach var of varlist ln_output {
		
		
			qui su `var' if hdmixd==0 & day_of_intervention==0
			local `var'refmean = r(mean)
			bootstrap, seed(17) cluster(line_team) strata(line)  nodots: ///
            reg `var' ldmixed* hdmixd? hdmixd?? hdmixd??? i.line i.Shift,  nocons 
			
			parmest, level(90, 95) fast 

			drop if strpos(parm,"ldmixed")==0 & strpos(parm,"hdmixd")==0
			g ldmixed = strpos(parm,"ldmixed")==1
			g day_of_intervention = substr(parm,8,.) if ldmixed==1
			replace day_of_intervention = substr(parm,7,.) if ldmixed==0
			destring day_of_intervention, replace
			replace day_of_intervention = day_of_intervention + 2 if ldmixed==0

			qui su estimate if ldmixed==1 & day_of_intervention==0
			local refdiff = r(mean)-``var'refmean'

			foreach var2 of varlist estimate min90 max90 max95 min95 {
				replace `var2'=`var2'-`refdiff' // shift so that mean of ldmixed on first bin is accurate
			}
			
			if "`var'"=="ln_output" {
				local yhead "Log (Output)"
			}
			
			
			qui su min90
			local yref = r(min)*1.01

					
* Refined Plot with Simplified Legend
twoway  ///
    (scatter estimate day_of_intervention if ldmixed==1, msymbol(circle_hollow) msize(medium) mcolor(orange)) /// LD-Mixed point estimates with hollow red circles
    (scatter estimate day_of_intervention if ldmixed==0, msymbol(square_hollow) msize(medium) mcolor(navy)) /// HD-Mixed point estimates with hollow blue squares
    (line estimate day_of_intervention if ldmixed==1, lcolor(orange) lwidth(medium) lpattern(solid)) /// Solid connecting line for LD-Mixed
    (line estimate day_of_intervention if ldmixed==0, lcolor(navy) lwidth(medium) lpattern(dash)) /// Dashed connecting line for HD-Mixed
    (rspike min90 max90 day_of_intervention if ldmixed==1, lcolor(orange%90) lwidth(medthick)) /// 90% CI for LD-Mixed with solid thin line
    (rspike min90 max90 day_of_intervention if ldmixed==0, lcolor(navy%90) lwidth(medthick)) /// 90% CI for HD-Mixed with solid thin line
    (rspike min95 max95 day_of_intervention if ldmixed==1, lcolor(orange%70) lwidth(thin)) /// 95% CI for LD-Mixed with very thin line, slightly lighter red
    (rspike min95 max95 day_of_intervention if ldmixed==0, lcolor(navy%70) lwidth(thin)), /// 95% CI for HD-Mixed with very thin line, slightly lighter blue
    title("`title'") ///
    xtitle("Days Since Intervention") ///
    ytitle("`yhead'") ///
    legend(order(1 "LD-Mixed Line" 2 "HD-Mixed Line") ///
           ring(0) pos(5) region(lstyle(none))) ///
    ylabel(9.9(0.2)11.1, nogrid) ///
    xlabel(0(20)100, nogrid) ///
    xlabel(0 "0-20" 20 "21-40" 40 "41-60" 60 "61-80" 80 "81-100" 100 "> 100", nogrid)
	
graph export "$Output/Figures/Event_Diag_AO1_BW.pdf", as(pdf) replace

					
	}
}



*--------------------------------------------------
* Robustness Tests (Permuation Test)
*--------------------------------------------------
if `ritest' == 1 {

    * Load Cleaned Data
    use "$Data/Final/linelevel_data.dta", clear
    sort line Date 
    
    /*--------------------------------------------------
    * Log Output in Pieces
    *--------------------------------------------------*/
    areg ln_output i.HD_Mixed i.line i.Shift, absorb(doy) cluster(line_team)
    estimates store modelA
    local obsA = e(N)
    su ln_output if e(sample) == 1
    local meanA = string(round(r(mean), 0.01), "%9.2f")
    
    * Perform RI Test for Model A
    ritest HD_Mixed _b[1.HD_Mixed], reps(2000) cluster(line_team) seed(9393) kdensityplot ///
        strata(line) nodots: areg ln_output i.HD_Mixed i.line i.Shift, absorb(doy) cluster(line_team)
    
    * Store p-values for Model A
    matrix pvaluesA = r(p)
    mat colnames pvaluesA = 1.HD_Mixed
    estimates restore modelA
    estadd matrix pvalues = pvaluesA
    
    /*--------------------------------------------------
    * Log Output in Boxes
    *--------------------------------------------------*/
    areg ln_productioncfc i.HD_Mixed i.line i.Shift, absorb(doy) cluster(line_team)
    estimates store modelB
    local obsB = e(N)
    su ln_productioncfc if e(sample) == 1
    local meanB = string(round(r(mean), 0.01), "%9.2f")
    
    * Perform RI Test for Model B
    ritest HD_Mixed _b[1.HD_Mixed], reps(2000) cluster(line_team) seed(9393) ///
        strata(line) nodots ///
        : areg ln_productioncfc i.HD_Mixed i.line i.Shift, absorb(doy) cluster(line_team)
    
    * Store p-values for Model B
    matrix pvaluesB = r(p)
    mat colnames pvaluesB = 1.HD_Mixed
    estimates restore modelB
    estadd matrix pvalues = pvaluesB
    
    * Export Robustness Test Results to LaTeX
    estout modelA modelB using "$Output/Tables/tables_line_level_ritest", style(tex) replace ///
        cells(b(star fmt(3)) p(par) pvalues(par([ ]))) ///
        nolabel collabels(none) mlabels(none) ///
        keep(1.HD_Mixed) ///
        starlevels(* 0.10 ** 0.05 *** 0.01) ///
        varlabels(1.HD_Mixed "HD-Mixed Line (LD Non-Mixed)")
    	
    * Create LaTeX Footnotes for RI Test Table
    local tex " \\ \hline"
    local tex "`tex' Day F.E. & Yes & Yes \\"
    local tex "`tex' Shift F.E. & Yes & Yes \\"
    local tex "`tex' Line F.E. & Yes & Yes \\"
    local tex "`tex' Mean Dep Var. & `meanA' & `meanB' \\"
    local tex "`tex' N & `obsA' & `obsB' \\"
    local tex "`tex' \multicolumn{3}{p{10cm}}{\tiny \textit{Notes:} Standard errors clustered at line-section-team level. P-values are reported in parentheses. Mixed is a dummy variable coded 1 if the line-section-level team is religiously mixed. *** p<0.01, ** p<0.05, * p<0.1.} \\ \end{tabular} }"
    
    * Append Footnotes to the LaTeX RI Test Table
    esttab modelA modelB using "$Output/Tables/tables_line_level_ritesta.tex", style(tex) replace booktabs ///
        d(*) nolabel noobs postfoot("`tex'") nonum ///
        mtitles("(1)" "(2)") ///
        mgroups("Log Output (Pieces)" "Log Output (Boxes)", pattern(1 1))
		
		
}



*--------------------------------------------------------
* Dynamics using binsreg command in STATA (older version)
*-------------------------------------------------------
if `dynamics_binsreg' == 1 {

    * Load Cleaned Data
    use "$Data/Final/linelevel_data.dta", clear
    sort line Date 

    * Create Bins for Days Since Intervention
    xtile bin_int_days = int_days, n(6)
    replace bin_int_days = bin_int_days + 0.15 if HD_Mixed == 1

    * Label Line Types
    generate line_type = ""
    replace line_type = "HD-Mixed Line" if HD_Mixed == 1
    replace line_type = "LD-Mixed Line" if HD_Mixed == 0
    
    /*--------------------------------------------------
    * Plot Dynamics Bins
    *--------------------------------------------------*/
    binsreg ln_output bin_int_days, absorb(line Shift) ///
        line(1 1) ///
        bylpatterns(dash) ///
        by(line_type) ///
        vce(cluster line_team) ///
        ci(1 1) ///
        level(90) ///
        bycolors(blue red) ///
        ytitle("Log Output") ///
        xtitle("Days Since Intervention") ///
        xlabel(1 "0-20" 2 "21-40" 3 "41-60" 4 "61-80" 5 "81-100" 6 "> 100") 
}
